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We present a simple implementation of the dynamical mean-field theory approach to the elec- 
tronic structure of strongly correlated materials. This implementation achieves full self-consistency 
over the charge density, taking into account correlation-induced changes to the total charge den- 
sity and effective Kohn-Sham Hamiltonian. A linear muffin-tin orbital basis-set is used, and the 
charge density is computed from moments of the many body momentum-distribution matrix. The 
calculation of the total energy is also considered, with a proper treatment of high-frequency tails of 
the Green's function and self-energy. The method is illustrated on two materials with well-localized 
4/ electrons, insulating cerium sesquioxide Ce203 and the 7-phase of metallic cerium, using the 
Hubbard-I approximation to the dynamical mean-field self-energy. The momentum-integrated spec- 
tral function and momentum-resolved dispersion of the Hubbard bands are calculated, as well as the 
volume-dependence of the total energy. We show that full self-consistency over the charge density, 
taking into account its modification by strong correlations, can be important for the computation 
of both thermodynamical and spectral properties, particularly in the case of the oxide material. 



2 



I. INTRODUCTION 

While density functional theory (DFT)^-^ in conjunction with the local density approximation (LDA) is 
remarkably successful in predicting ground-state properties of a wide range of real materials, it has been 
found unable to provide the correct description of so-called strongly correlated materials (transition metal 
oxides, many actinide and lanthanide-based materials, high Tc superconductors) even on a qualitative level. 
In order to overcome these shortcomings of the traditional DFT-LDA scheme, a combination of the DFT- 
based band structure techniques with dynamical mean-field theory (DMFT}^ has been proposed^i^i^. In 
DMFT one introduces a strong local Coulomb interaction acting between electrons of a correlated band (for 
example, the d— band in transition- metal oxides or the /—band in actinides). The self-energy is found by 
first mapping a full solid onto a quantum impurity model involving a single atom hybridized with an effective 
bath, and solving this effective model using many-body techniques. This self-energy is then promoted to all 
lattice sites, therefore restoring the translational invariance of the crystal. As a result one obtains the fully 
interacting Green's function for the system, from which a wide range of properties can be extracted. 

Starting from the end of the 90 's, this new LDA-I-DMFT approach to the electronic structure of strongly 
correlated materials has been rapidly developing, and a number of different implementations have been 
proposed and applied to calculations of the spectral and, in some cases, thermodynamic properties, of Mott 
insulators, ferromagnetic 3-d metals, Ce, Pu, and other actinide systems (for reviews, see^i^i^,). Up to date 
most LDA-I-DMFT calculations have been performed using a partially self-consistent scheme, where the 
local self-energy is obtained from a DMFT calculation with the fixed LDA charge density, and, hence, with 
a fixed LDA Hamiltonian. Therefore, in this simplified scheme one neglects the impact of the strong on-site 
Coulomb interaction on the charge distribution. 

Implementing full self-consistency over the charge density is a somewhat delicate task. Several fully 
self-consistent LDA-I-DMFT schemes have been discussed in the literaturei^iiiii^ii^, while only two actual 
implementations have appeared up to date. The first one, due to Savrasov and Kotliar— , is based on the 
full-potential LMTO method. In order to compute the LDA-I-DMFT charge density one has to construct 
from the LDA Hamiltonian and the DMFT local self-energy a local Green's function of the interacting 
system. In Ref. ^ this task has been accomplished by first finding the (right- and left-) eigenfunctions 
and eigenvectors of the Kohn-Sham Hamiltonian combined with the self-energy HksO^) + S(iw) and then 
performing the Brillouin zone integration by means of the tetrahedron method^^. As the self-energy is a 
complex function, this scheme requires the diagonalization of a non-Hermitian matrix for each k— point and 
Matsubara frequency, which is a computationally demanding task. Another scheme has been proposed by 
Minar et alM- on the basis of the KKR Green's function method. Their approach requires solving the radial 
Dirac equation with the self-energy being added to the LDA atomic potential. Generally the self-energy is 
a non-diagonal matrix, therefore an additional coupling is generated between Dirac equations for different 
magnetic quantum numbers m, making the solution of the Dirac equation a highly non-trivial task. Minar 
et al. applied their scheme to calculations of the FeNi binary alloy, which forms a cubic lattice. Hence, 
in this case, both the Green's function and self-energy are diagonal in the orbital indexes, and there is no 
additional coupling generated in the Dirac equation. However, this problem will certainly arise for lower 
symmetry lattices or /—electron compounds. 

In the present article, we propose and apply a relatively simple implementation of the fully self-consistent 
LDA-I-DMFT method. The charge density (including correlation effects described by the many-body self- 
energy) is expressed in terms of a k-dependent momentum-distribution matrix N^j^,, which is obtained from 
the many-body Green's function by summing over frequency. Furthermore, using a linear muffin-tin orbital 
basis set, we reduce the calculation of the charge density to that of three moments, involving this matrix and 
the Kohn-Sham Hamiltonian itself. In order to calculate the total energy, a functional of both the charge 
density and the on-site components of the Green's function associated with local orbitals can be used^ii^, 
as previously discussed by Savrasov and Kotliar. In all these calculations, obtaining the charge density and 
total energy with sufficient accuracy requires a careful treatment of the high-frequency tails of the Green's 
function and self-energy, and we derive here appropriate formulas to handle this issue. 

In order to illustrate this approach, and to assess the importance of full self-consistency on the charge den- 
sity, we perform in this article fuUy-selfconsistent LDA-I-DMFT calculations of the density of states, the band 
structure, and the volume-dependence of the total energy of two materials: 7-cerium and cerium sesquioxide 
Ce203. The quasi-localized /-electrons arc treated using DMFT in conjunction with the Hubbard-I approx- 
imation as a quantum impurity solver^^. Our theory reproduces the experimentally observed splitting of the 
rare-earth /-band into occupied lower and empty upper Hubbard bands, while the conventional LDA cal- 
culations incorrectly predict this band to be pinned at the Fermi leveU^. We show that the self-consistency 
over the charge density shifts significantly the positions of the Hubbard bands (in comparison to using a 
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frozen LDA charge density), in both 7-Ce and Ce203. It is therefore important for the correct description 
of spectral and optical properties of those materials. 

The paper is organized as follows. In Section HH we describe our self-consistent LDA+DMFT implemen- 
tation. After a brief reminder of the basic LDA-I-DMFT scheme (Sec. Ill Ap and of the LMTO basis-set 
fSec. IIIBj) . we derive the expression for the charge density in this basis fSec. Ill C|) . After reviewing the total 
energy functional of LDA+DMFT (Sec. Ill D|) . we discuss the practical calculation of the energy (Sec. IIIE|) . 
We conclude Sec. |lT]by describing the specific form of the interaction vertex used in this article, as well as 
the double counting correction and the Hubbard-I impurity solver (Sec. lIIF"]) . Section IIIIl presents the results 
for CezOg (Sec. and 7-Ce (Sec.|niB|). 

II. IMPLEMENTATION OF THE FULLY SELF-CONSISTENT LDA+DMFT METHOD 
A. The LDA+DMFT formalism: a brief reminder 

The LDA+DMFT approach to electronic structure is based on two key quantities: the charge density p(r) 
and the local Green's function of the solid, or more precisely the projection Gati^^) of the full Green's function 
onto a single atomic site and on the subspace of correlated orbitals. Both quantities have to be determined 
self-consistently, following an iterative cycle which is summarized on Fig. [1] This can be rationalized as a 
functional of both p(r) and Gab, as detailed later in Sec. Ill Dl 

Let us follow this iterative cycle, starting from a charge density profile p{r). A Kohn-Sham (KS) potential 
is constructed from p{r) as: 

VKsi'r) = Vc{r) + vh [p{r)] + -w^^jc [p(r)] (1) 

in which Vc is the periodic potential of the (fixed) ions, vh = J dr' j^^^p(r') is the Hartree potential, 
and Vxc = SExc/Sp{r) is the exchange-correlation potential. Hence, the functional dependence of vks on 
p(r) is kept identical to that of conventional DFT. In practice, Vxc will be computed from p(r) using the 
LDA form of the exchange-correlation energy E^^"^ = J drp{r)e[p{r)], with e[p] the energy density of the 
homogeneous electron gas. Solving the single-particle Schrodinger equation associated with vks yields the 
KS eigenenergies e^^f and eigenfunctions |kt^) (with v a band index), forming the KS effective one-particle 
hamiltonian: 

Hks = 5]6^,^|ki.)(kH (2) 

At this stage, it is convenient to introduce a set of localized basis functions XLR(r), where R denotes an 
atomic position, and L stands for all orbital indices (e.g., L = {I, m, a}). In the following, we shall consider 
linearized muffin-tin orbitals (LMTOs), but different basis-sets can be used and have been considered in the 
literaturei^iiilii^ (e.g., Wannier functions). The electron creation operator at a point r in the solid can be 
expanded on this basis as: 

^t(r)=ExlRW4R (3) 

R,L 

and the KS hamiltonian reads: 

HkS ^Y.^L3i^)clLC^L' (4) 
kL 

In the LDA+DMFT approach, this Hamiltonian is supplemented by many-body terms. These many- 
body terms act in the subspace generated by a set of orbitals corresponding in practice to the orbitals (e.g 
d or /-orbitals) for which a description beyond DFT-LDA is needed (note however that all other orbital 
components in the valence will also be modified indirectly by feedback effects of the self-energy associated 
with the correlated ones). For simplicity, we restrict the discussion here to the case of one 'correlated' 
atom per unit cell. The orbitals generating the correlated subset need not coincide with basis functions in 
generali^. However in the present work, we do choose them as a specific subset XaR(i") of the LMTOs. Hence, 
L = {1,771, a} runs over all orbitals retained in the valence, while a = {m, cr} runs only over the 'correlated' 
subset (denoted by C in the following). The many-body Hamiltonian considered in LDA+DMFT reads: 



H = Hks — Hdc + Hu 



(5) 
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In this expression, Hoc (corresponding to a one-body potential Vdc) is a double-counting correction. Indeed, 
some of the local Coulomb interaction effects are already taken into account in the exchange-correlation 
energy, and hence in Hks- The many-body terms Hjj act in the subset of correlated orbitals only. They 
correspond to matrix elements of the Coulomb interaction, and will in general involve general 2-particle 
terms of the form ^^^^ U abcdC^a'RAiH.'^d'R.Cc-R.- (In the present work however, only density-density terms 
will be retained: the form of Hjj and Hdc used in this article is discussed further in Sec. Ill F( ). 

Let us consider the full Green's function of the solid G(r, r'; r — t') = — (TV'(r, t)%P'^{y\ t')), which can be 
decomposed on the basis set as: 

G(r, r'; r^T') = J2J2 ^^i^(^) ^^^'^^ - R', r - r') XL'R'(r')* (6) 

RR' LL' 

DMFT focuses on the local components Gab of the Green's function, on the same atomic site (R = R') and 
within the correlated subspace. The key idea (which can be viewed as a representability assumption) is that 
Gab can be represented by an effective local model, which is a multi-band generalization of an Anderson 
impurity model described by the effective action: 



S^ - [\t dr' J2 [Go'Uir - t') c,(r') + I dr Hu (7) 

ab 



"'0 



In this expression, Qq is the dynamical mean-field, analogous to the familiar Weiss mean-field in classical 
mean-field theory, the key difference being that here it is a frequency-dependent (i.e energy-scale dependent) 
quantity. It can also be viewed as the hybridization function which connects the correlated atom (effective 
impurity) to its environment: Aab{z) — zSab ~ — [Qo{z)^^]ab- In this expression, z is an arbitrary 
frequency in the complex plane, and e/ is a matrix of effective on-site atomic levels (see sec. IIIFI) . The 
dynamical mean-field Qq (or A) is determined from a self-consistency condition, which expresses that the 
impurity-model Green's function faithfully represents the local Green's function in the solid projected onto 
the correlated subset, and hence that the two quantities should coincide. Furthermore, an approximation is 
made, namely that the many-body self-energy has components on the basis set which are (i) local and (ii) 
non-zero only in the correlated subspace, so that it also coincides with its impurity model counterpart and 
takes the form: 



^LL^ (2) = f^R.R' ( n ) (^) 





s:7(.) 



The impurity model Green's function and self-energy are defined as: 

Gl7(r - r') ^ -(Tct (r)c,(r'))„„p , K7 ^ [Q^'U - [Gi^.U (9) 

in which the average indicated by {■■■)imp is taken with respect to the effective action ([7]). The self-consistency 
condition which determines Qq can thus be expressed in a concise way as: 

-Fr^-Pr = Gimp (10) 

in which the full Green's function G of the solid is projected on a given correlated atom, and onto the 
correlated subset with the projector: 

a 

Given ([9]) and given Dyson's equation relating G and E, Eq. pO)) yields an implicit relation between Gimp 
and Qq. To be fully explicit, let us go through the iterative process (the DMFT loop depicted in Fig.[T|) which 
determine Gimp^'^imp and self-consistently. Given an initial guess for Q^^ the impurity Green's function 
Gimp is calculated using some appropriate 'impurity solver', and the impurity self-energy is obtained as: 
E^™^ = [t/(7^]a6 — [G~mj^ab, and used as the only non-zero block of T^ll' (Eq. [5]). The Green's function in 
the solid can then be calculated as: 

[G-i]ii,(k, z) = (z + ^l)SLL' - H^3 + - ^ll\z) (12) 
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and summed over k-points in the Brillouin zone in order to yield the local Green's function on the atomic 
site associated with the correlated atom as: 

Giociz) = E [(^ + - + VES ^Lv{z)\ (13) 

k 

The block corresponding to the correlated orbitals is extracted from this expression and inverted in order to 
get an iterated value for the dynamical mean-field, as: 

[g-Xb = s„„p + [p^GiocP^]-' , (14) 

and the DMFT iteration can be pursued. Note that, even though the self-energy matrix has only components 
in the subspace of correlated orbitals, components of the Green's function corresponding to all valence orbitals 
(s,p, • • • ) are modified due to the matrix inversion. Hence there is a feedback of correlation effects onto these 
orbitals as well. 

At the end of the DMFT cycle, the resulting Green's function can be considered as a full many-body 
solution of the Hamiltonian defined by Eq. ([5]), for a given charge density determining H^g, under the 
approximation that the corresponding self-energy is taken to be purely local. The charge density can then 
be recalculated from the Green's function at the end of this cycle, as: 

= E E ^^i^W - R', r - r' = 0-) XL'R'(r)* (15) 

RR' LL' 

From this updated charge density, a new Kohn-Sham potential is constructed from ([T]), and the associated 
one-particle Schrodinger equation is solved again in order to produce an updated H^g, which serves as 
a new input for the DMFT cycle. This is indicated as the DFT loop on Fig. [T] Full self-consistency is 
reached when all local quantities {G imp, imp, Go) as well as the charge density p(r), have converged. Note 
that the LDA-I-DMFT self-consistent charge density is affected by correlation effects through the many-body 
self-energy and that it differs from its LDA self-consistent value plda{j^)- Furthermore, it also differs from 
the charge density associated with the occupied Kohn-Sham orbitals and evaluated with the converged Hks 
but E = in (jisp. i.e: pks{y) = X^ki/ IV'kiy(r)|^. That is, the Kohn-Sham representation is not used for the 
self-consistent charge density within LDA+DMFT. In the next two sections, we first give a brief reminder 
of the LMTO formalism and basis-set, and then proceed with the practical evaluation of expression ^5]) for 
the charge density in this basis-set. 



B. The LMTO basis set 

As discussed in detail in section III A( since DMFT emphasizes local correlations, one needs to construct a 
localized basis set, i.e basis functions which are centered on the atomic positions R in the crystal lattice. Up 
to now, most implementations have used basis sets based on linear muffin-tin orbitals^^ (LMTOs) XiR(i") = 
XL(r — R). These basis sets offer the advantage to carry over the physical intuition of atomic orbitals from the 
isolated atom to the solid. The LMTO method has been extensively used in electronic structure calculations 
and thoroughly described in review articles2£i2i. Here we only outline the main features of the LMTO basis 
set which are relevant to the selfconsistent implementation of the LDA-I-DMFT scheme. 

MT-methods suppose that the crystal potential is spherically symmetric near each atomic site and constant 
in the interstitial region between atoms. The conventional LMTO method employs a further simplification, 
the atomic-sphere approximation (ASA), which assumes that the whole space of a crystal is filled with 
atomic spheres and neglects both the overlap and interstitial spaces. The LMTO basis is constructed from 
the solution (partial wave) 0lr of the spherically symmetric KS potential inside the atomic sphere located 
at R for a certain energy E^, (typically, the center of gravity of a band) and its energy derivative 0lr- 
The angular dependence is provided by a spherical harmonic Yj, , corresponding to the orbital and magnetic 
quantum numbers L). The expression for a linearized MT-orbital is: 

X2R(r) = 0LR(r) + E (r)/i2'R'LR, (16) 

L'R' 

where the first and second terms in the right-hand side of this equation are usually called the "head" and 
" tail" , respectively. The superscript a designates a particular LMTO representation (the so-called "screened" 
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FIG. 1: (Color online). DMFT combined with electronic structure calculations. Starting from a local electronic 
density p(r), the associated Kohn-Sham potential is calculated and the Kohn-Sham equations are solved. The Kohn- 
Sham hamiltonian Hj^f,{]i) is expressed in a localised basis set (e.g LMTOs). A double-counting term is subtracted 
to obtain the eflFective one-electron Hamiltonian Hq = H''"^ — H^^ . The local self-energy matrix for the subset 
of correlated orbitals is obtained through the iteration of the DMFT loop: a multi-orbital impurity model for the 
correlated subset is solved (red (lighter grey) arrow) , containing as an input the dynamical mean-field (or Weiss field 
Qo). The self-energy Tiah is combined with Hq into the self-consistency condition Eqs. (|13I14|I in order to update 
the Weiss field (blue (darker grey) arrow). At the end of the DMFT loop, the components of the full, k-dependent, 
Green's function in the local basis set can be calculated, yielding the momentum-distribution matrix 7V(k)£,j^/ and 
the updated charge density p(r) as described in Sec. Ill CI This updated charge density is used to compute the new 
Kohn-Sham potential until a converged local density is also reached (for the DFT loop). At the same time, the 
chemical potential must also be adjusted self-consistently with the total number of electrons. 



non-orthogonal representation) , which is defined by the choice of the envelope functions for the interstitial. 
We have actually employed here the so-called nearly-orthogonal 7-representation2^, where is chosen as: 

Using the orthogonality property of the partial wave and its energy derivative ((/>lr|(/)lr) = one may easily 
show^^ that the head of the LMTO in the 7-representation is orthogonal to any LMTO centered on any other 
site up to a second-order contribution due to the overlap of the energy derivatives of the partial waves. Then 
the basis can be made completely orthogonal by means of an numerical orthogonalization, for example, 
the Lowdin transformation^. The use of an orthogonal basis is technically simpler when implementing 
LDA-I-DMFT (although a non-orthogonal basis may be more advisable in principle since it is expected to 
reduce the range of interaction terms). A numerical orthogonalization usually introduces undesirable mixing 
between strongly and weakly correlated states. Because the LMTO basis in the 7 representation is already 
nearly-orthogonal, the inter-orbital mixing are small however. 

One of the main advantages of the LMTO technique is the small size of its basis. It can be made even 
smaller through use of the downfolding procedure^!, which allows to reduce the size of the Hamiltonian 
by folding down those states which are located well above the valence band. This technique substantially 
reduces the computational effort, a gain which is especially important in the case of the time-consuming 
LDA-HDMFT calculations. The down-folded LMTO can be expressed as follows: 

xIr(i") = 0LR(r) + 0L'R'(r)^I'R'LR + 0ifR" (i-)^ifR"LR, (18) 
Z,'R' h-r." 

where L labels "active" orbitals, which are explicitly included in the Hamiltonian, H runs over the set 
of downfolded orbitals, the matrix Z can be expressed through the structure constants and potential 
parameters^!. Following Ref. Q we neglect the energy dependence of the downfolded orbitals cj^HR" (r) • In 
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the case of a periodic crystal, one may rewrite expression (fT8|) at each k-point in the Brillouin zone as: 

X^(r) = Mr) + E ^L'{r)h^L'L + E '^^W^^l. (19) 

L' H 

C. The Calculation of the Charge Density 

In this section, we describe the practical implementation of the calculation of the LDA+DMFT charge 
density in the LMTO basis set. Expression (fT5|) can be rewritten as: 

/^W =EE^^k(r)A^^^,xI'k(r) (20) 

LL' k 

in which N^j^, is the k-dependent occupancy matrix related to the full Green's function ^ by: 

7V^^, = GLL.(k,r = 0-) =T^GLL'(k,ic^„)e'"°^ (21) 

n 

The last equality is expressed as a sum over fermionic Matsubara frequencies cj„ = (2n + l)7r//? corresponding 
to the temperature T — 1/(3. By inserting the linear MT-orbitals (fT9|) in ([20|). one can obtain a general 
formula for the LDA+DMFT charge density in the LMTO framework, which involves simply the calculation 
of momentum averages of N and of products of N with the hamiltonian matrix h. This expression can be 
further simplified if the atomic-sphere approximation (ASA) is used, in which the potential is approximated as 
spherically symmetric within the MT-spheres. In accordance with this approximation, crossed terms between 
different angular momentum channels can be neglected. Neglecting also the overlap between spheres and 
using the orthogonahty = and {4'l\<I>h) = between partial waves yields our final expression 

for the angular-averaged charge density in a given unit cell (valid within the ASA): 

p{r) = J2 «^\Mr)\' + 2mfc^L{r)k{r) + mf\Mr)\^) + J] ('') (22) 

L H 

In this expression, the moments to^*-''s are defined from the k-dependent occupancy matrix and LMTO 
first-order Hamiltonian h (for a given atom R) as: 

mf ^5]iV^^ = (trN)k 

k 

<^ = E E Z^hlNIl'ZIh = (tr(ZHNZH))k (23) 

k LL' 

m«^^^iV^^,/i|,^ = (tr(Nh))k 

k L' 

^T.T. hlvNhv'hln^ = (tr(hNh))k 

k L'L" 

Expressions (|22|23|1 are similar to those used in the context of usual DFT in the LMTO-ASA formalism, the 
key difference being that in the DMFT context, the momentum-distribution matrix N'^ is computed from 
the many-body Green's function according to (|20[) instead than from filling independent orbitals as in the 
KS representation of the density. In practice, the moments ([23|) are computed at the end of the DMFT 
cycle and then passed on into the LMTO electronic structure part of the program, where a new total charge 
density is computed according to formula as indicated on Fig. [T] 

D. The Total Energy Functional 

In order to discuss total energy calculations in the LDA-I-DMFT framework, it is best to use a formu- 
lation of this scheme in terms of a (free-) energy functional. Kotliar and Savrasov^SiS^i^ have introduced 
for this purpose a ("spectral-density-") functional of both the total charge density p(r) and the on-site 
Green's function in the correlated subset: (denoted Gab for simplicity in the following). Let us 
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emphasize that these quantities are independent, since Gab is restricted to local components and to a 
subset of orbitals so that p(r) cannot be reconstructed from it. The functional is constructed by intro- 
ducing source terms Z(r) — vks{^) ~ Vc{r) and A'Sab{i^^n) coupling to the operators tp'' {r)ij;{r) and to 
^j^X*(r — R)i/'(r, r)'0'''(r', r')x6(r' — R) — Ca'R{T)clj^{T'), respectively. Furthermore, the Luttinger-Ward^i 
part of the functional is approximated by that of the on-site local many-body Hamiltonian Hjj — Hue 
introduced above. This yields: 



Sl[p(r), Gab; VKs(y)-, AS a\^LDA+DMFT = 

-tr ln[zcj„ + /i + - VKs{^) - X*-AS .x] - J dr {vks - Wc)p(r) - tr [CAE] 
-Hi / dr dr'p(r) ^ p(r') + E^Mr)] + Er i^.mp[GT] - ^dc[GT]) 



+ 



In this expression, x*.AS denotes the "upfolding" of the local quantity AS to the whole solid: x*-AI] = 
Sr Safe — R-)Sa6(«a;„)xf)(r' — R)- Variations of this functional with respect to the sources Jfi/J = 
and 6il/6Y,ab — yield the standard expression of the local density and local Green's function in terms of 
the full Green's function of the solid, which we have used in the previous section: 



p{r) = (r|G|r} , Gab{iuJn) = {Xa-R\G\xbR) 



with: 



G = 



2^'-^'A-s(r)-X*-AE.x 



(24) 



(25) 



Note that these expressions, as well as the functional, are written here in a manner which does not refer 
explicitly to a specific basis set. The formalism does depend, however, on the choice of localized orbitals 
defining the correlated subspace. 

From these relations, the Legendre multiplier functions vks and AE can be eliminated in terms of p and 
Gab, so that a functional of these local observables only is obtained: 



^lda+dmft[p, Gab] = ^LDA+DMFT [p(r), Gafc; l[p, G], AE [p, G]] (26) 

Extremalisation of this functional with respect to p {6T /5p — 0) and Gab {6^/5 Gab = 0) yields the expression 
of the Kohn-Sham potential and self-energy correction at self-consistency: 



"2 5E. 



.K5(r)=..(r)+/dr'^— ^p(r') + ^ (27) 



^^ab-^Ti =^ab - ^ab 

OLrab OUab 



Hence, one recovers from this functional the defining equations of the LDA+DMFT combined scheme detailed 
in the previous section, including self-consistency over the local density ([M)) . 
The free-energy (|24l26p leads to the following expression of the total energy: 



Elda+dmft = Ek,LL' H^3(y)NLU +Idr[v,{r) - VKs{rMr) + (29) 
+ i / drdr'pir) ^ p(r') + E,,,[p] + {Hu) ~ Edc- 

The sum of the first four terms in this expression is reminiscent of the expression of the total energy in DFT, 
Edft[p{^)], evaluated at the self-consistent charge density p(r), except for the fact that the many-body 
(LDA-I-DMFT) momentum-distribution matrix enters the first term. Instead, within DFT, the first term 
reads X^ki/ ^kif ^ where the prime indicates that the sum is to be taken only over occupied KS orbitals (filled 
up to pks, adjusted so that the total number of electrons is obtained. Hence, in practice, the total energy 
can be calculated as: 

Elda+dmft = ( i?Di=^T[p(r)] - ^ etA + J2 H^HmL' + (Hu) - Edc = 

V ku / k,LL' 

= E, [p] + Eh [p] + S.c [p] + i^LL' + {Hu)- Edc (30) 

k,LL' 
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The first three terms are the crystal, Hartree and exchange-correlation energy, respectively. The last three 
terms coincide with the energy associated with the many-body Hamiltonian Hxs—Hjyc+Hjj. the interaction 
energy {Hjj) may be computed either directly from the expression of the Hamiltonian Hu (for example, by 
evaluating the correlations (nafih) within the impurity solver, which is easy e.g when using quantum Monte 
Carlo). Alternatively, the Migdal formula (Hu) — Tr {YjG)/2 can be applied, as actually done in this work 
when using the Hubbard-I approximation and described in the next section. 



E. Some practical aspects of total energy calculations 

An accurate evaluation of the k-dependent occupancy matrix (|2ip is necessary in order to obtain the 
correct charge density (|22|) and total energy ([29]) . If calculations are performed at finite temperature then 
this requires a careful summation of the high-frequency tails of the Green's function in Ij2ip . In addition, 
the formula (|29p for the total energy contains the Migdal contribution, which is also computed through the 
corresponding summation over the Matsubara frequencies: 

{Hu) - ^Tr [g(t = 0-)S(r = 0")] - ^Tr ^ G(zc^„)I](zu;„)e^'^"0^ (31) 

n 

Here both the Green's function and self-energy contain high-frequency tails, which should be properly taken 
into account in the evaluation of the sum over LOn- In Refs. [28ll29ll30t evaluation of the Matsubara sums 
over the high-frequency tails have been treated in details, however the approach proposed there is applicable 
only if the self-energy is calculated by means of an analytic technique. Due to the wide spread of numerical 
DMFT solvers (for example, quantum Monte Carlo), it is highly desirable to have a "solver-independent" 
technique for accurate evaluations of Matsubara sums, which we describe in this chapter. 

We shall isolate the first two terms in the high-frequency expansion of the self-energy matrix, which is 
thus decomposed as: 

A ~ 

Y.{iLOn) = i;(ioo) + h T.num{iUJn) = T,an{iuJn) + S,mm(«W„) (32) 

In this expression, the numerically determined S„„m(za'„) = S(ia;„) — Ea„(ici;„) will be neglected for Mat- 
subara frequencies larger than a certain cutoff LUcut, while the high-frequency contribution Sa„(iw„) = 
I](ioo) -|- A/iujn will be treated analytically when performing frequency sums. One may extract the matrices 
S(joo) and A (the latter being actually diagonal in the present case) from the real and imaginary parts of 
the self-energy at the cutoff frequency E(iWcMt), or use a more sophisticated way for fitting them. A similar 
decomposition can be applied to the Green's function: 

G(k, iUn) = G„„„i(k, iuj„) + GaniK iuJn) (33) 



G„„™(k, iujn) + [iuJn + - ^ + y - S(ioo 



-1 



where Graum(k, iun) = G(k, iuun) — Gan(k, iujn) is again zero for Matsubara frequencies larger than the cutoff 
frequency. In the analytical part of the Green's function it is sufficient to keep only the dominant term 
in the self-energy E(ioo). Then the matrix /i -I- V^'~^ — H^^ — E(ioo) is Hermitian, and we designate its 
eigenvectors and eigenvalues as \X^-^) and Aj^, respectively, so that: 

G(k,zw„) = G„„™(k,ic^„)-f V^^^^. (34) 

In order to evaluate the frequency sum in (|2ip we carry out the summation of Gnum up to the cutoff 
frequency, while other frequency sums are performed analytically. By inserting (j32p and (|34p in the formula 
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for the Migdal energy (|3T|) one obtains the following expressions of practical use: 

(Hu) = (Hu)'^'^ + {Hu)'^'^ + {Hu)'^^\ (35) 
T r 

num {iojn) + G num 

(k, za;„)I](iw„) 

k ni 
k m 

The Hubbard-I quantum impurity solver— employed in the present work is not constrained to the Mat- 
subara frequencies, and it can in fact be equally well used for calculating the self-energy at any general 
complex energy. Hence at zero temperature one may easily rewrite frequency sums in the expressions for the 
occupancy matrix (|2ip and Migdal energy (|29p in a form, which is suitable for a summation over the poles 
of the GF (and the self-energy in the case of Migdal energy) on the real axis: 

N^L' = ^jGLL'iKz), (36) 



(Hu) ^^^fTr [G(z)I](z)] dz, (37) 

where integration is performed over the contour in the complex energy plane, which encloses valence-band 
energy poles. Therefore, within the Hubbard-I approach and at zero temperature one may completely avoid 
the problem of summation of the high-frequency tails. We have applied both the finite temperature Mat- 
subara temperature summation and contour integration techniques in the fully self-consistent LDA-I-DMFT 
calculations of 7-Ce and Ce203. 



F. Choice of interaction vertex, impurity solver, and double-counting 

The last term in Eq.([5]), Hu gives the many-body interaction terms acting in the subset of correlated 
orbitals. They correspond to matrix elements of the Coulomb interaction, and will in general involve arbitrary 
2-particle terms Uabcdc\c\cdCc. There the 4-index matrix Uabcd is defined for a /-shell by four Slatter integrals 
F^, F^, and . In the quasiatomic (spherical) approximation the Slatter integrals can be expressed 
through only two parameters U and J (see Ref. i31i). In addition to the spherical approximation one often 
makes a further simplification and keeps only density-density interactions (we employ this simplification in 
the current version of our Hubbard-I impurity solver, though in general it is not necessary). We shall limit 
ourselves here to this case and use: 

^u = \T. E Kun'n^nAn (38) 

R, rnm'aa' 

with the effective two-index matrix derived from a more general 4-index form as follows: 

^ mm' — Ur/rfim'mm' ; ^ mm' — ^ mm' — U^am'mm' Ujji^' rn' ^ (39) 

Next, we discuss the impurity solver that we use in practice in this article. Both materials that we 
shall consider for illustrative purposes (Ce203 and 7-Ce) have nominally an configuration, and the /- 
electron is actually localized. The practical solution of the DMFT equation can be simplified considerably 
by choosing an approximate 'impurity solver' appropriate to this localized character. The simplest of those 
is the 'Hubbard-I' approximation. In this approximation, the self-energy is approximated by its 'atomic 
limit', in which the hybridization function Aab{z) is neglected. One should still correctly identify however 
the effective atomic levels entering the Weiss dynamical mean-field ['y(7"'^]mm' — z — emm' ~ Amm'{z) (with 
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m,m' G C). For this purpose, we can perform a high-frequency expansion of the self-consistency equation 
(|13ll4p and request that A(z) vanishes at high frequency, which leads to: 

ernm' = ^mrn' + (k) " K^m' (40) 

k 

In the present case, this is actually a diagonal matrix of effective atomic levels Emm' = ^m^mm'- These levels 
must be recalculated iteratively, as they change upon updating the chemical potential. Hence, the effective 
atomic Hamiltonian reads: 

= + (41) 

This Hamiltonian can be diagonalized, yielding (many-body) energy levels Ea's and eigenstates \A), from 
which the atomic Green's function can be constructed as: 

The atomic self-energy is obtained from S^t — (z — em)Smm' — G'^f . This leads to the following expression 
for the Green's function of the full solid, in the Hubbard-I approximation: 

[G(k,z)-i]ii, = + + V''')5lu - H^3i^) + [G-/(z)]ll. (43) 
^[G-,\z)]LL'-H^3ik) (44) 

in which H^^, = Hff,{k) - {Hff,{k))k for L,L' eC and H^f, ^ H^3ik) otherwise (i.e. H^'^ is the non- 
local, inter-atomic part of the KS Hamiltonian). Within the Hubbard-I approximation, the LDA-I-DMFT 
loop takes the following form. Starting from H^^ at a given stage of the iteration, the effective atomic 
levels are calculated according to (|40|) and the local atomic Green's function calculated by diagonalizing the 
effective atomic Hamiltonian (HI]). The full Green's function is then formed as and the momentum 
distribution matrix N^j^, calculated. The chemical potential is then updated so that ^ N^j^ yields the 
appropriate total number of electrons, and the updated charge density is obtained by calculating the moments 
as described in Sec. Ill CI The KS equations are then solved to yield a new H^^ and this process is iterated 
until convergence of both p{r) and of the effective atomic levels Cm- 

Finally, we discuss the "double-counting" correction Hoc- This correction must be introduced, since the 
contribution of interactions between the correlated orbitals to the total energy is already partially included 
in the exchange-correlation potential derived from E^c- Unfortunately, it is not possible to derive this term 
explicitly, since the energy within DFT is a functional of the total electron density, which combines all 
orbitals in a non-linear manner. In practice, the most commonly used form of the double-counting term is 
(for other choices, see e.g [s^): 

Hdc = E K.t^'cLcfc. 



*^ aba 



U{Nf-\)~J{N]^\) 



(45) 



where Nf = Nj + is the total occupancy in the correlated shell C (i.e the /-shell in practice in this 
article). One issue arises here, which is which value of Nf must actually be used in (|45|) . When solving 
the DMFT equations with a numerically exact solver, it would seem (from the previous functional-based 
derivation) that Nf should be the occupancy of the /-shell obtained at self-consistency. However, we have 
found this to be inappropriate when using the Hubbard-I approximation as a solver and leading to too small 
equilibrium volumes. Instead, we note that the Hubbard-I solver treats an effective isolated atom, which has 
this a frozen occupancy taking integer values at T = 0. Hence, a natural choice within Hubbard-I is to choose 
this frozen integer occupancy in the double-counting correction. For the /^ materials in the paramagnetic 
phase treated in this article {nJ — Nj — 1/2), so that the Hund's coupling therefore drops out of Hjjc 
which reads simply: 

U 
2" 

This also implies that the [/-dependent contribution of the double counting correction in the total energy 
actually vanishes for compounds within the Hubbard-I approximation, and for this choice of double- 
counting potential. 



12 



III. FULLY SELF-CONSISTENT LDA+DMFT CALCULATIONS OF CE2O3 AND 7-CERIUM 

We have applied the LMTO-based fuUy self-consistent LDA+DMFT technique in order to calculate the 
density of states and thermodynamical properties of the cerium sesquioxide Ce2 03 and pure Ce in its 7- 
phase. The main aim of these calculations is to validate our implementation of the fully self-consistent 
LDA+DMFT technique, as well as to study the impact of the charge self-consistency on spectral and ther- 
modynamic properties. Hence, it is desirable to avoid additional complications due to the use of sophisticated 
and computationally expensive quantum impurity solvers. While the Ce 4/ band cannot be properly mod- 
eled by conventional LDA techniques, due to its essentially atomic-like (localized) character, the electronic 
correlations on the 4/ shell can be treated by means of the simplest "strong-coupling" Hubbard-I (HI)^^ ap- 
proximate impurity solver, described in the previous section. Therefore for our purposes these two strongly 
correlated rare-earth compounds, with nominally one localized /-electron, appear to be a suitable choice. 
Both pure cerium and the oxide Ce203 are well studied theoretically and experimentally, because they raise 
questions of fundamental interest (e.g. concerning the 01-7 transition in 0^^'^*^ ) as well as for technological 
reasons (for example, Ce-oxides are used for oxygen storage in solid-oxide fuel cellsi^). 

A. CezOg 

The electronic properties of Ce203 are largely determined by the localized Ce 4/ orbitals. This compound 
is an insulator with a gap of about 2.5 eV. At low temperatures, Ce203 orders antiferromagnetically, with 
a Neel temperature (9 K) which is several orders of magnitude smaller than the gap. It is thus clear that 
the antiferromagnetic order is not the driving force behind the insulating character of this material. Rather, 
strong correlations open up a gap among the 4/ states, and this compound can therefore be viewed as 
an f-electron based Mott insulator. As for all such materials, it is a challenge for conventional electronic 
structure calculations to describe the opening of the correlation-induced gap. The Kohn-Sham spectrum of 
DFT-LDA is metallic, with 4/ states at the Fermi level. An insulating state of Ce203 has been obtained 
within the self-interaction correction approachS^, albeit with a strongly overestimated value for the band 
gap. In the antiferromagnetic phase, the LDA+U method can be employe d^^'^^i'^^'^^ , as well as the hybrid 
functional approach^. In this article, we focus on the paramagnetic phase, at low temperature just above 
the Neel ordering temperature. Dynamical mean-field theory and the Hubbard-I approximation used here 
allows us to describe the opening of the Mott gap and the existence of local moments in the paramagnetic 
phase, due to the quasi-localized 4/ electrons. 

1. Details of the calculational setup 

We have calculated the density of states and spectral density of Ce203 oxide in the hexagonal lattice 
structure (space group P3ml) at experimental volume (a =3.890 A, c/a— 1.557— ). Empty spheres have 
been introduced in order to make the structure more close packed. In the computation of the spectra we 
have included the 6s, 5d and 4/ orbitals of Ce, the 2p orbitals of oxygen and the Is orbitals on the empty 
spheres in the basis set as active LMTOs. In addition, the 6p orbital of Ce, the 3s and 3d oxygen orbitals, and 
the 2p orbitals on the empty spheres have been downfolded. In the total energy calculations it is necessary 
to include the Ce semicore 5p orbitali^ as well. The LMTO code employed by us is not able to treat semicore 
states in a separate panel, hence in the total energy calculations we have excluded the Ce 6p orbitals from 
the basis. It is important to include the 6p orbitals in order to obtain accurate spectra; however, as will be 
shown in the next chapter on 7-Ce, for total energy calculations they are less significant. 

In the self-consistent LDA+DMFT calculations we have employed 62 k— points in order to carry out the 
integration over the irreducible Brillouin zone of the hexagonal lattice, the energy integration has been 
carried out on a semicircular contour of depth 2 Ry comprising 40 energy points. In order to obtain the 
local and k-resolved spectral functions, the self-energy and Green's function were directly computed on the 
real axis, with 394 k— points used for the integration over the irreducible Brillouin zone. 

In order to estimate the value of the screened on-site Coulomb interaction we have performed constrained 
LSDA calculations^ both, by fixing the /—electron occupancy and by adding a constrained potential acting 
on the /-electrons on one Ce atom. Depending on whether the constrained charge or potential approach 
was used we obtained quite different values of the parameter U, varying from 5.5 to 8 eV. The previous 
theoretical estimates of Refs. of the U value in pure Ce metal put it in the range between 5 and 6 

eV, while in the LDA+U calculations of Refs. (sslfsgj the best agreement between the calculated equilibrium 
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FIG. 2: (Color online). The LMTO Ce2 03 DOS calculated within DFT-LDA, the LDA+DMFT with the fixed 
LDA charge density (DMFT-nonSC), and the LDA+DMFT fully self-consistent over the charge density (DMFT-SC) 
schemes. The black (solid), green (dashed), blue (dotted) and red (dash-dotted) curves are the total, O, Ce— d and 
Ce— / DOS, respectively. The vertical dashed line indicates the position of the chemical potential, while the vertical 
dotted lines indicate positions of the lower Hubbard bands in DMFT-nonSC and DMFT-SC, respectively. 



volume and experiment was achieved for a value of the effective Ucff — U — J of about 6 eV. Thus, in our 
calculations we have fixed the value of [/ — J to 6 eV, while the value of J=0.46 eV, which depends less 
crucially than U on the proper treatment of screening in a material, has been taken from the constrained 
LDA calculations. 

We have employed three calculational schemes, for comparison purposes: the conventional LDA, the 
(non self-consistent) LDA-hDMFT with fixed LDA charge density (DMFT-nonSC), and the LDA+DMFT 
technique with full self-consistency over the charge density (DMFT-SC). In the case of DMFT-nonSC, we 
carried out DMFT iterations until convergence in the self-energy and chemical potential is reached. In the 
case of DMFT-SC, additional convergence criteria with respect to the total energy and charge density were 
employed as well. 

2. Spectral functions and correlated 'bandstructure ' 



The orbitally-resolved local density of states (DOS) - or spectral fmictions-, defined as: 

Al{uj) = --ImV GLL(k,c^ + iO+) (46) 

77 ^ — ' 
k 



14 



KS 



{Hi) Vdc jej) 



DMFT-nonSC 



1.19 



1.52 



1.18 -3.00 -3.01 



DMFT-SC 



1.15 



2.76 



2.18 



-3.00 



-1.97 



TABLE L The chemical potential, parameter C of the KS / band, double counting and / level positions in the 
DMFT-nonSC and DMFT-SC approaches (in eV). 

are displayed in Fig.[2]for all three methods (LDA, DMFT-nonSC, and DMFT-SC). First, one may notice that 
conventional LDA calculations place the Ce 4/-band in the vicinity of the Fermi level, therefore predicting 
Ce203 to be a metal. Both the DMFT-nonSC and DMFT-SC approaches correctly predict Ce203 to be a 
Mott insulator, where the /-band is split due to the local Coulomb interaction into occupied lower Hubbard 
bands, and empty upper Hubbard bands. One may notice a significant shift of the positions of the Hubbard 
bands in the DMFT-SC DOS with respect to the DMFT-nonSC picture. The value of the band gap in 
Ce203 is equal to 3.10 eV and 2.13 eV within the DMFT-nonSC and DMFT-SC approaches, respectively 
The latter value is in good agreement with the experimental measurements of the optical gap in Ce203 
(2.4 eV, Ref. [3]), while the fixed charge (non-SC) calculations lead to a strong overestimation of the 
gap. The total occupancy of the /-shell is rather weakly affected by the local Coulomb interaction. The 
occupancy of the / shell in Ce203 is equal to 1.145, 1.174, and 1.167 according to the LDA, DMFT-nonSC 
and DMFT-SC calculations, respectively. 

In order to better understand the observed difference between the position of the lower Hubbard band in 
the non-SC and SC methods, we use the expression (^0]) of the effective atomic levels within the Hubbard-I 
approximation, established in the previous section. This yields directly the position of the lower Hubbard 
band (and of the upper Hubbard band, which is shifted upwards by an energy U) since, for both materials, 
the / shell of the isolated atom contains just one electron. The effective level position for the m-th orbital 
component of the 4/ shell reads: 



The double counting term used in the present work reads: Vdc — J7(iV/ — 1/2) — J(iV/ — l)/2 and since we set 
Nf = 1 in this expression (as explained in the previous section), we see that the double counting correction 
term is the same for both the DMFT-nonSC and DMFT-SC approaches. Thus, the difference in the value 
of the band gap cannot be blamed on the double counting. Therefore, in the Hubbard-I approximation, and 
for a given local Coulomb interaction, the factor which determines the positions of the lower Hubbard band 
is the momentum average of Hks, i-e essentially the centre of gravity of the KS f-band (the second term in 
the right-hand side of expression (|T7|) ). 

In Table U we list the chemical potentials, the LMTO parameter of the Ce / band (which can be 
interpreted as the center of a "pure" (unhybridized) KS band), the Hamiltonian integrated over the BZ, 
the value of the double counting correction and the level position. The brackets (...) in Table U designate 
an average over the magnetic quantum number m within the f-shell. It is important to realize, that while 
in the case of DMFT-nonSC, the parameters C^g and are just obtained within the conventional LDA, 
in DMFT-SC it is obtained from the solution of the KS equations with a potential calculated from the 
self-consistent LDA-I-DMFT charge density, which differs from the LDA one. One may notice that the KS 
band structure computed from the DMFT charge density is substantially different from the LDA one, as 
the centerweight of the / band is shifted by 1.24 eV upwards. This change in the KS band structure is 
reflected in the values of the Hamiltonian (H^), and, finally, in the positions of the DMFT Hubbard bands, 
which are determined by the level position e,„. Hence one observes a clear impact of the charge-density 
self-consistency on the LDA-I-DMFT electronic structure: changes in the charge density due to DMFT lead 
to a modification of the KS band structure, which in turn affects the final LDA-I-DMFT electronic structure. 
The result is a reduction of the band gap by almost 1 eV. A direct connection between the KS and DMFT 
electronic structures is especially evident in the present case due to the simplicity of the HI technique, where 
the position of the LHB is linked to the KS Hamiltonian through the simple relation (|T7)) . 

Since the Hubbard-I approximation employed here does not describe lifetime effects of the correlated states, 
the momentum-resolved Green's function still has well-defined poles (in other words, excitations correspond- 
ing to the lower and upper Hubbard bands have infinite lifetime) . The dispersion of the corresponding upper 
and lower Hubbard bands of Ce203 as a function of momentum is displayed in Fig. [31 as obtained with fully 
self-consistent LDA+DMFT. The LHB shows virtually no dispersion and forms a narrow, atomic-like level 
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FIG. 3: (Color online). The LMTO fully self-consistent LDA+DMFT k-resolved spectral function of Ce203. 

in the gap between the oxygen p and Ce d band. The UHB is essentially a set of atomic- like "multiplet" 
states, while one may still notice some evidence of weak hybridization between the d band and the UHB^^ . 

3. Total energy calculations 

The total energy calculations have been carried out by keeping the experimental c/a ratio and using 
the same setup as in the calculations of the spectra, apart from modifications in the basis, where the 5p 
Ce semicore states were included instead of the 6p orbitals. Without the 5p semicore states there is no 
minimum of the total energy vs. volume curve in the range of lattice parameters between 3.60 and 4.2 A. 
As in the calculations of spectra, we have employed both the fully self-consistent and fixed LDA charge 
density implementations of the LDA-I-DMFT method. In order to compare our results obtained within the 
LMTO-ASA with the recent full-potential calculations of Andersson et al^, we have also carried out LDA 
and LDA-fU calculations of Ce203. The LDA-fU calculations have been performed using the same setup 
and U and J parameters than for LDA-I-DMFT, but using a Hartree-Fock approximation to the self-energy 
in the self-consistent calculations. The HF self-energy is obtained for an / occupancy, which can be different 
from the atomic one, therefore, in contrast to the HI case, in the LDA-I-U calculations we have used the 
actual (self-consistent) /-band occupancy in the expression for the double-counting correction. 

Our results for the energy vs. lattice parameter curves are displayed in Fig [3) One may see that the 
localization of the Ce 4/ electrons due to the strong local Coulomb correlations leads to a substantial 
increase in the lattice parameter. This is expected on a physical basis, since localization leads to a decreasing 
participation of these electrons to the cohesion of the solid. Our results for the equilibrium lattice parameter 
within the DMFT-nonSC, DMFT-SC and LDA+U methods are rather close to each other (3.79, 3.81 and 
3.84 A), and in reasonable agreement with the value 3.86 A, obtained in Ref. [3^ for a similar value of U, as 
well as with experiment (3.89 A). The LDA lattice constant is much smaller (3.72 A). The difference between 
the DMFT-nonSC and DMFT-SC for the equilibrium volume appears, somewhat surprisingly, to be rather 
small. The difference between the LDA-fU and LDA-I-DMFT equilibrium volumes (i.e between Hartree-Fock 
and Hubbard-I approximations for the self-energy) stems possibly from the non-conserving nature of the HI 
approximation. 



16 




FIG. 4; (Color online). The total energy vs. lattice parameter dependence calculated within LDA (solid line), 
LDA+U (dashed line), DMFT-nonSC (dotted line) and DMFT-SC (dash-dotted line) approaches.The vertical dashed 
line indicates the experimental lattice parameter of Ce2 03 



B. 7-Cerium 




CO (eV) CO (eV) 

FIG. 5: (Color online). Spectral function for 7-Cerium within LDA+DMFT (Hubbard I) at T=OK with the self- 
consistent (SC) and non self-consistent (nonSC) schemes. In the first case, two ways of treating the double counting 
corrections are used: LDA number of electron and DMFT one (N=l for the Hubbard I solver). Experimental 
spectrum (PES and BIS) of 7 Ceriuro^i^ are also represented. 



As a second example, we have applied our DMFT-SC code to fcc-Cerium. This system has been thor- 
oughly studied in the past in order to understand the a-^ transition'^'^, an isostructural volume collapse 
transition taking place as a function of temperature. Above the transition temperature Tc, in the large 
volume (7-) phase, the f-electrons are localized, while the smaller volume in the a-phase below Tc leads to 
a more delocalized behavior of the f-electrons. Several model studies have been carried out to describe this 
transition^i^i^. More recently, it has also been studied in the framework of ab-initio calculations using 
LDA-l-DMFTi^i^i^iSi^i^i^i^. However, these studies use the fixed charge density (non-SC) implemen- 
tation of LDA-I-DMFT. Nevertheless, they were successful in describing the main aspects of the transition 
(spectra, and Kondo stabilization energy). The goal of our study is to evaluate the importance of self- 
consistency effects in LDA-f DMFT calculations for Cerium. Since the most recent calculation o^^i^^'^^ were 
done with the precise but time consuming Quantum Monte Carlo, it was not possible to test it. Another 
important goal is to study if a simple Hubbard 1 implementation of LDA-I-DMFT is able to describe correctly 
the 7 phase of Cerium, in which electrons are more localized than in a Cerium. 

Computations were done at the experimental volume of the 7 phase for the spectra (34. SA ). 5p semicore 
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FIG. 6: (Color online). Total energy versus lattice parameter curves in LDA+DMFT (Hubbard I) for Cerium (shifts 
in energy between different curves are arbitrary). 

states were taken into account thanks to an implementation of the LDA+DMFT self-consistency over density 
scheme with the multiple LMTO code^^. Calculations with 5s states included do not show any change in 
the results. Valence states contain 5p 6s 6p 5d and 4f states. Calculations were done with 145 k-points in 
the irreducible Brillouin zone. 

Figures [H] and [5] shows the spectral function and the total energy versus volume curves for fcc-Cerium 
for both the non-SC and the SC calculations. In accordance to calculations on Ce203, we use the 'atomic' 
number of electron computed in the Hubbard I solver to evaluate the double counting correction. It is more 
justified because of error compensation and, as we show below, the results for the lattice parameter at T=0 
are consistent with LDA+U self-consistent calculations using the self-consistent number of electrons in the 
solid. For testing purposes, a non-SC calculations has also been done with the LDA number of electron in 
the double-counting correction. 

Concerning the spectral function, we first see that all the calculations give the same qualitative physical 
pictures. The distance in energy between Hubbard bands is the same in all the calculations. The difference 
lies in the exact position of the Hubbard bands. The SC calculation gives the best agreement with experiment, 
the non-SC is shifted by 0.2 eV. The non-SC calculation with LDA double-counting is shifted by 1 eV with 
respect to the SC calculation. Table [H] gives the data necessary to understand these shifts, as explained 
before. We see that the effect of self-consistency itself is small in this case. The main error in the non- 
SC/DClda comes from the double counting term, as expected. Note that non-SC and SC calculation lead 
to nearly the same number of electron. This can be partly attributed to the correlated character of the 7 
phase of Cerium: electrons are localized, so the number of electron is close to 1 (at least compared to an 
LDA calculation) . Experimental spectra^^i^ of the gamma phase are also represented on figure O While 
the calculated position of the Hubbard bands is in good agreement with experiment, the width of the upper 
Hubbard band is not correctly described in our calculations. This is not surprising and is to be blamed on 
the use of a simplified Hubbard-I solver, which docs not include lifetime effects. This means that Hubbard-I 
is only partly adequate to describe the high temperature 7-phase. 

Figure [6] gives the total energy versus volume computed at T=OK. The main conclusion is that non-SC and 
SC DMFT calculations give the same minimum for the internal energy versus volume curves (see also table 
mil) . This value is not far from the value obtained in the LDA+U calculations in ASA (which is itself in good 
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DMFT-nonSC/DCLDA 


1.06 


0.23 


0.45 


3.94 


-3.72 


DMFT-nonSC/DCN=i 


1.04 


0.26 


0.45 


2.99 


-2.80 


DMFT-SC/DCn=i 


1.04 


0.10 


0.57 


2.99 


-2.52 



TABLE II: The number of electrons in DMFT, the chemical potential, the double counting contribution to the 
potential and the / level position in the DMFT-nonSC and DMFT-SC approaches for Cerium. Energies are in eV. 
The number of f-electrons in LDA is 1.16. 
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agreement with LDA+U calculations in PAW^-) . Note however that LDA+U can only account for correlation 
effects by introducing a spurious magnetic order, while DMFT is able to describe local moment formation 
and correlation effects in the paramagnetic phase. Effect of self-consistency appears to be fairly weak in the 
case of 7-cerium. One should also keep in mind that the ASA approximation does not describe correctly 
the bulk modulus in LDA. So the results concerning the energy versus volume curves should be taken only 
as trends. Moreover, direct comparison with experiment is difficult, because we have not computed the full 
free energy but only the internal energy^. 



Expffi£l 5.16 
LDA+U/PAW^ 5.05 
LDA+U/ASA 5.00 
DMFT SC/ASA 4.93 
DMFT non-SC/ASA 4.91 
LDA/PAW^ 4.52 
LDA/ASA 4.52 

TABLE III: Lattice parameter (in A) of 7 Cerium according to experimental data, PAW calculations and our calcu- 
lations. 

The calculations where 6p states are not taken into account leads to an underestimation of the lattice 
parameter (4.80 A instead of 4.93 A). This is in accordance to the underestimation of the lattice parameter 
in LDA when 6p states are neglected (4.44 A instead of 4.52 A) 

Our conclusion regarding 7-cerium is that spectra and energy are reasonably described with a fixed (LDA, 
non-SC) charge density approximation. However, it is mandatory to use a self-consistent implementation to 
obtain more accurate results. 



IV. CONCLUSION 

In conclusion, we have devised a simple and efhcient implementation of the fully self-consistent 
LDA-I-DMFT method in the LMTO basis set. The charge density is calculated from moments involving 
the LDA-I-DMFT momentum-distribution matrix and the KS hamiltonian. We have also obtained accu- 
rate formulas for computing the total energy by handling high-frequency tails of the Green's function and 
self-energy in an appropriate manner. 

We have computed the local and k-resolved spectral functions of cerium sesquioxide Ce203 by means of 
the fully self-consistent LDA-I-DMFT technique in conjunction with the Hubbard-I approximation for the 
DMFT self-energy. We have shown that the charge-density self-consistency affects the spectral properties of 
Ce203 substantially, causing a shift of the lower Hubbard band by approximately 1 eV and a corresponding 
decrease in the value of the band gap in comparison with DMFT calculations with fixed LDA charge density; 
This effect considerably improves the agreement with experiment in comparison to a non-self consistent 
calculation using the LDA charge density. We have identified the main cause of these modifications, which is 
due to a significant change in the effective KS band structure, computed with LDA-I-DMFT charge density, 
as compared to the LDA band structure. 

Finally, we have obtained the total energy and equilibrium volume of 7-Ce within the fixed (LDA) charge 
and fully self-consistent LDA-I-DMFT schemes. The effects of the self-consistency over the charge density 
are less important for 7-Ce than for the oxide: self-consistency induces in this case a change of 1% on the 
ground-state volume, and affects the spectral function by shifts of a fraction of an electron- Volt only. 
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